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We analyze the main features of granular shear flow through experimental measurements in a 
Couette geometry and a comparison to a locally Newtonian, continuum model of granular flow. The 
model is based on earlier hydrodynamic models, adjusted to take into account the experimentally 
observed coupling between fluctuations in particle motion and mean-flow properties. Experimentally, 
the local velocity fluctuations are found to vary as a power of the local velocity gradient. This can be 
explained by an effective viscosity that diverges more rapidly as the random-close-packing density 
is approached than is predicted by Enskog theory for dense hard sphere systems. Experiment and 
theory are in good agreement, especially for the following key features of granular flow: The flow 
is confined to a small shear band, fluctuations decay approximately exponentially away from the 
sheared wall, and the shear stress is approximately independent of shear rate. The functional forms 
of the velocity and fluctuation profiles predicted by the model agree with the experimental results. 
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I. INTRODUCTION 

The general features of granular shear flow have been 
investigated thoroughly over the last several years [§],§||. 
The following key points emerge in shear flow experi- 
ments on a large range of materials in two and three 
dimensions: 

• The velocity of particles decreases quickly over a 
few particle diameters away from the shearing wall 
(see e.g. §-§). 

• The velocity profile, normalized by the shear rate 
U, is independent of U (see e.g. f|||). 

• The shear force a is approximately independent of 
U, if the granular material is allowed to dilate (see 
e.g. 1). 

These features, together with the discovery of strong 
inhomogeneities in the force distribution even during 
flow , might be taken to indicate that any continuum 
approach, such as local hydrodynamic models, should fail 
to describe granular flow. 

Here we revisit the assumptions made in earlier hydro- 
dynamic models of granular flow, via careful comparison 
to the experimentally measured microscopic particle dy- 
namics in a circular Couette geometry. This leads us to 
emphasize the strong interplay between local RMS fluc- 
tuations, the mean flow, and the local density. When this 
coupling is properly taken into account, a hydrodynamic 
model, which we have introduced in Ref. H, quantita- 
tively describes all key properties of granular shear flow 
discussed above, including both flow properties and shear 
forces. 

The shear force obtained from the hydrodynamic 
model resembles the simple dynamic friction law found in 



solid-on-solid friction, i.e. the shear force is proportional 
to pressure, but approximately independent of shear rate. 
We emphasize that this result is obtained even though the 
hydrodynamic model does not include frictional forces 
between grains. 

The temperature is defined here as the mass times the 
square of RMS velocity fluctuations. While the temper- 
ature is a constant in shear flow of an ordinary fluid, 
granular temperature is dissipated through inelastic col- 
lisions. Its spatial variation plays a crucial role in de- 
termining the properties of granular flow. We find that 
the granular temperature profile, normalized by its max- 
imum value, is independent of shear rate and pressure. 
Temperature is introduced into the system via viscous 
heating over a characteristic length of the order of a few 
particle diameters. It is then dissipated via inelastic col- 
lisions over a longer length scale. Our model contains 
a description of the source and transport of fluctuations 
that allows us to predict the pressure and the shear rate 
dependence of both the shear forces and the particle dy- 
namics. 

In the experiments reported here, which go beyond 
those reported earlier [Q , the granular material is sheared 
in a Couette geometry with a rotating inner cylinder and 
a stationary outer cylinder. The inner cylinder is con- 
nected to the motor through a flexible spring that allows 
either stick-slip motion or steady shearing depending on 
parameters. We can also apply an upward air flow at a 
variable rate through the granular material to dilate the 
material and reduce the stresses. Our ability to vary the 
stresses via air flow, and our study of both stick-slip and 
steady dynamics in the same apparatus, distinguish this 
work from other experimental measurements of sheared 
granular matter H,||. 

In addition to performing force measurements, which 
probe macroscopic material properties of the ensemble of 
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particles, we determine the dynamics of individual par- 
ticles by measuring the mean velocity and RMS velocity 
fluctuations of particles on the top surface of the granular 
layer. The combination of velocity and force measure- 
ments, together with variation of the stresses and the 
time dependence of the flow, allow for a very sensitive 
test of our theoretical model. 

Previous experimental results and modeling ap- 
proaches are reviewed in section [ij. Our experimental 
setup and results for the particle dynamics and shear 
forces in sheared granular matter are discussed in sec- 
tions [II and IV. In section [y], the locally Newtonian 
hydrodynamic model is described in detail. We conclude 
in section VI with a discussion of the main results and 



the broader implications of this work. 



II. BACKGROUND 

Efforts to understand particle dynamics during granu- 
lar flow can be roughly divided into adaptations of con- 
tinuum models, based on hydrodynamic or elasto-plastic 
descriptions, and models that emphasize the differences 
between molecular fluids or solids and granular matter, 
such as the inhomogeneous character of particle contacts 
and of stress transmission in a granular material. In this 
section we review briefly these different approaches. We 
conclude with a discussion of the aim of our hydrody- 
namic model in the context of previous work. 



A. Continuum models of granular flow 

Hydrodynamic models were motivated by Bagnolds 
pioneering theoretical and experimental work on shear 
forces in dense suspensions The constitutive equa- 
tions for macroscopic quantities, such as the shear stress 
as a function of shear rate, were investigated in several 
studies of dense suspensions with fixed volume. In the 
limit of large velocities, Bagnold found that in a fixed 
volume the shear stress a is proportional to the shear 
rate U squared. He referred to this regime as the "grain 
inertia regime" . He accounted for the measurements by 
assuming that the local shear stress a was proportional 
to the square of the local shear rate 7. He justified such a 
relationship on the basis of kinetic arguments. Since di- 
rect measurements of microscopic particle dynamics were 
not available, he assumed a linear velocity profile similar 
to that of ordinary fluids. 

The development of hydrodynamic descriptions, ex- 
tracted from the "microscopic" dynamics using kinetic 
theory, was pioneered by Jenkins and Savage || and by 
Haft pc[ |. A considerable amount of theoretical, numeri- 
cal, and experimental work has refined and modified this 
approach. Reviews of models of granular flow that de- 
scribe many of these studies has been compiled by Camp- 
bell jnj, Savage and more recently by Clement H. 



In the kinetic theory approach, the granular material 
is generally modeled as an inelastic hard-sphere system. 
Constitutive equations similar to the usual Navier-Stokes 
equations of hydrodynamics can be obtained in the limit 
of small inelasticity. The transport coefficients entering 
the flow equations are usually computed at the level of 
the Enskog equation J35|] , an extension of the Boltzmann 
equation that takes the finite size of particles into account 
but neglects correlations between collisions. Due to the 
assumptions that enter into the kinetic theories presented 
above, these descriptions are limited to rapid granular 
flows and to intermediate or low particle densities. 

In order to understand the boundary between a seem- 
ingly flowing state and an apparently stationary state, 
various extensions of the previous hydrodynamic model 
have been proposed. Jenkins and Askari 1 13 have studied 
the interface between a flowing region and an amorphous 
(high density) region which is at rest. In this model, the 
thickness of the shear band is determined by the balance 
between the energy input and the loss rate due to inelas- 
ticity On the other hand, visco-plastic models have been 
proposed. The Savage-Hutter model |l4| uses the Mohr- 
Coulomb failure criterion to predict the transition from 
solid-like to fluid-like behavior in the context of avalanch- 
ing and rock slides. The constitutive relation connecting 
shear stress to shear rate is of great practical importance. 
Different relations have been proposed for various situa- 
tions, and several are summarized in a table in [jl5| . 



B. Alternative descriptions of granular flow 

Beyond measurements of the mean properties (such as 
shear forces) that are important for a continuum model, 
recent experimental and theoretical studies have focused 
on measurements of the particle dynamics and shear 
forces at the scale of an individual particle. 

Detailed measurements of the particle dynamics (see 
e.g. [§-[|) revealed that in several experimental geome- 
tries particle motion is confined to several (generally 5- 
10) particle diameters close to the sheared surface. The 
velocity profile is found to be roughly exponential or 
Gaussian. 

The force distribution within a granular assembly, 
measured with birefringent disks p6| or carbon pa- 
per |l7j ], was found to exhibit strong inhomogeneities on 
the particle scale. Stresses were found to be transmitted 
along chains of particles (force chains) in a static granular 
assembly and during shear. 

There have been several attempts to account specifi- 
cally for these inhomogeneities in granular flows. Some 
approaches describe the flow properties on the basis of 
fracture models jl8|;[l9| , while others introduced non-local 
constitutive equations coupling force chains to flowing 
grains f§. 

The shear strength of deformable, inelastic spheres was 
modeled using a discrete element method by Aharonov 
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and Sparks j6| . The density is found to adjust to a critical 
density within the shear band. 

C. Our hydrodynamic model 

While theories and models of granular flow have be- 
come remarkably detailed, the assumptions of each model 
can strongly influence the results. As noted by Camp- 
bell [ fi"l| , detailed experimental measurements of particle 
velocities, granular temperatures and densities were often 
unavailable when models were developed. The measure- 
ments reported here should help assess the validity of the 
assumptions of various models. 

Here we revisit the local hydrodynamic model and 
carefully reexamine the assumptions made to derive the 
constitutive equations in view of our experimental mea- 
surements of individual particle dynamics and mean 
shear forces. 

The mean shear forces have been studied before in an 
experimental system similar to ours: Tardos et at. §] 
investigated the effect of an upward air flow through the 
granular material on the shear forces. The shear force 
was found to decrease linearly with air flow. On the basis 
of Bagnold's results || and much subsequent work, it is 
clear that dilatancy has an important effect on granular 
flow. This was also convincingly demonstrated in the 
experiments by Tardos et ai: If the material is allowed 
to expand, the shear stress is independent of shear rate. 
On the other hand, if the material is confined to a fixed 
volume, Bagnold's result of a quadratic increase of shear 
stress with shear rate is found. 

The interplay between dilatancy, fluctuations, and flow 
is complex. Theoretical investigations have not yet pro- 
duced results that are consistent with measured mean 
(macroscopic) properties and measured dynamics at the 
particle (microscopic) scale. Granular material develops 
a greater resistance to flow as its density increases. Pre- 
vious theoretical treatments of granular flow have either 
been restricted to the lower-density rapid-flow regime 
jjj, or they have incorporated a yield threshold (visco- 
plastic models) that produces a well-defined transition 
from liquid- like to solid-like behavior [Q . 

Recent experimental studies of granular flow down a 
sandpile by Nasuno et al. [pi] using long exposure time 
video imaging have revealed that the transition between 
solid-like and fluid-like behavior may not be very well 
defined. The velocity profile within the flowing layer was 
found to be exponential over more than seven orders of 
magnitude in speed with no clear transition to a solid-like 
state. These results indicate that no strict transition may 
exist, but instead one may be able to treat the solid- like 
state as a very high viscosity fluid. 

In our model we do not distinguish between a solid-like 
and fluid-like state. Rather, we try to model increased 
resistance to flow in the high-density limit within the 
standard kinetic approach (i.e., we assume binary colli- 



sions, no friction, etc.). Previous approaches pjl0| used 
the Enskog equations to take excluded volume at higher 
densities into account. However, correlated motion of 
particles, e.g. through cooperative rearrangements, is not 
taken into account in the Enskog equation. 

The high density regime of polydisperse spheres is 
bounded by the largest possible random close packing 
(RCP), which experimentally is found to be 63.7% for 
slightly polydisperse systems p^ |. For a wider distri- 
bution of particle sizes, slightly larger densities can be 
reached. In order to model the high density regime 
close to RCP, we use the results from elastic systems 
as a first approximation. In elastic systems, simulations 
in Lennard-Jones systems [^3| and experiments on col- 
loids |^6| suggest that cooperative rearrangements are 
important at high density. In these systems the viscos- 
ity is found to diverge anomalously strongly with density 
in the high density limit p4jj25|. The divergence with 
density is roughly equivalent to the divergence with de- 
creasing temperature for a supercooled liquid close to the 
glass transition |2(| . 

The introduction of a viscosity that diverges more 
rapidly than the thermal conductivity and heat loss co- 
efficient as random close packing is approached is the 
key new feature of our model. It leads to theoretical 
predictions in quantitative agreement with most of our 
experimental results. It also provides a direct connection 
between the dynamics of granular media and glasses, as 
has been proposed by Liu and Nagel [^7) . 

III. EXPERIMENTAL SETUP 
A. Apparatus 

In the experiments we shear the granular material 
in a Couette geometry. The granular material used 
in most of the experiments reported here consists of 
0.55 — 0.95 mm diameter black glass beads (from Jaygo 
Inc.) (p m — 2.55 g/mm 3 ). The color does not alter sur- 
face properties, but increases the opacity of the material, 
which facilitates the tracking of particles on the surface 
as described below. We also carried out experiments with 
a mixture of 1.3 mm and 1.6 mm Chrome Steel spheres 
(Pm = 5.0g/mm 3 ), and with polydisperse, rough ceramic 
spheres (Macrolite ML1430 from Kinetico Corp.) with 
diameter 0.83 — 1.47 mm (p m = 0.51g/mm 3 ). 

In the experimental apparatus the granular material 
is confined to a 12 mm gap between a stationary outer 
cylinder and a rotating inner cylinder (r = 51 mm), as 
shown in Fig. [3. The gap can be reduced to 3 mm. The 
inner cylinder is hollow to reduce its inertia and is coated 
with a monolayer of randomly packed glass beads to pro- 
vide a rough boundary. The outer glass cylinder is coated 
with a monolayer of randomly packed glass beads up to 
the height of the top surface, which allows observation 
of the top layer of grains through a mirror as shown in 
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Fig. |. The lower 38 mm of the inner cylinder is station- 
ary in order to minimize boundary layer effects. 

To shear the material, the inner cylinder is ro- 
tated with a 4000 step/turn microstepping motor (from 
Aerotech Inc.) at a variable rate of 0.001 — 1 Hz. The 
rotation rate is smoothed by a 100 : 1 gearhead for rates 
< 1 Hz. The inner cylinder is connected to the microstep- 
ping motor via a flexible tempered steel spring. This 
spring configuration allows us to measure instantaneous 
shear forces with excellent dynamic range and precision, 
since the spring bending is proportional to the applied 
shear force. We measure the spring bending with a ca- 
pacitive displacement sensor (EMD1051, Electro Corp.) 
that is rigidly connected to the motor shaft at a radial 
distance of 4.2 mm from the shaft. The spring constant 
of the spring (dimensions: 0.51 x 7.5 x 165 mm) was de- 
termined to be 220 ± 8 N/m. 
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FIG. 1. Experimental setup: The granular material (be- 
tween two concentric cylinders) is fluidized by an upward air 
flow and sheared by rotation of the inner cylinder, which is 
connected to the motor through a flexible spring. Shear forces 
are determined from the spring displacement. Particle mo- 
tions in the top layer are measured through the glass outer 
cylinder with a fast CCD camera. 

The soft connection between the motor and the inner 
cylinder permits both stick-slip dynamics and continuous 
motion of the inner cylinder to be obtained, depending 
on parameters. However, when a uniform speed of the 
cylinder is required for the experiment, the spring is re- 
placed with a rigid connection. In that case measurement 
of shear forces is not possible. 

We can also apply an upward air flow at a variable 
rate through the granular material. The air flow enters 
the granular material through a circular opening between 
the cylinders, and leaves the granular material through 
a circularly symmetric opening at the top. This assures 
a uniform air flow rate throughout the material. Flow 
rate uniformity was tested by observing the position of 
air bubbles as they leave the upper surface of the granu- 
lar material. Air bubbles form at high flow rates in the 



class of granular materials used in this experiment. The 
random position of air bubbles indicates that the flow 
is uniformly distributed throughout the gap. Except for 
this test for flow uniformity, none of the experiments pre- 
sented in this paper was carried in the presence of air bub- 
bling. The air flow is provided by a blower (Rigid Inc.) 
operated at variable input AC voltage. Flow rates are 
measured by means of an air velocity transducer (FMA 
904, Omega Inc.). For the range of air flow rates we 
employed, the average density of the granular material 
changes by less than 10 % and the flow speed is calcu- 
lated in first approximation assuming the porosity for 
random close packing. 

Since the air flow exerts a drag on individual grains, 
the effective weight supported by neighboring grains de- 
creases with increasing upward airflow. This effective 
weight of individual grains in turn should be proportional 
to the pressure inside the granular material. We can, 
therefore, reduce the pressure by applying an upward air 
flow, and increase pressure by applying a downward air 
flow. The proportionality factor between air flow and 
pressure can be roughly estimated by calculating the up- 
ward drag exerted on a single sphere at the mean air flow 
speed within the granular material pa] . 



B. Determination of particle dynamics 

We measure the mean particle velocities V(y) and 
the velocity fluctuations 5V{y) on the upper surface of 
the granular material. These should approximate parti- 
cle motion in the interior based on a previous measure- 
ments || that found very similar velocity profiles in the 
interior (measured with MRI and X-ray techniques) and 
on the bottom surface of a shear cell. We return to this 
issue in Sec. VI. 

The trajectories of individual particles in the surface 
layer are determined with a fast CCD camera at 30 — 
1000 frames/sec. Particle motion is extracted from four 
sequences of 2000 images using procedures written in IDL 
(RSI Inc.) based on tracking routines provided by J. 
Crocker and E. Weeks. 

In the first step of the tracking process, both long 
range brightness fluctuations (i.e. nonuniform illumina- 
tion) and short wavelength noise are reduced by applying 
a bandpass filter with a short wavelength noise cutoff and 
a long wavelength cutoff of roughly one particle diameter. 
In the second step, the position of particles (roughly 100 
in each image) is determined by calculating the centroid 
of each bright region in the filtered image. This yields 
a spatial resolution of < 0.1 pixels, provided that the 
bright region is several particle diameters wide. Black 
glass beads are better suited than undyed glass beads 
for an accurate determination of particle positions, since 
black beads are more opaque. This reduces internal re- 
flections and reflections from particles in deeper layers. 
In order to improve spatial resolution, the intensity peaks 
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are broadened by taking images slightly out of focus. The 
broader intensity peaks improve the precision of the cen- 
troid determination. Defocusing also reduces the inten- 
sity of secondary peaks due to scattering by multiple par- 
ticles to a level where they are no longer interpreted as 
particles. For the ceramic particles the defocusing pro- 
cess eliminates multiple peaks due to the substructure of 
individual particles. 

In the third step of the trajectory determination, the 
particles are labelled and the evolution of their position 
through an image sequence is determined. The assign- 
ment of particles to corresponding points in the previous 
and next frame is based on a tracking algorithm, which 
minimizes the total squared displacement within a se- 
quence of frames. In a final step, the probability dis- 
tribution of individual particle displacements is used to 
verify that large displacement particles are not system- 
atically cut off. 

Since the mean-flow velocity V is comparable to the 
RMS velocity fluctuations close to the inner cylinder, ac- 
curate tracking of particles is only possible if the maxi- 
mum displacement is considerably smaller than the dis- 
tance between particles. We have verified that particles 
are accurately tracked, even when the maximum parti- 
cle displacement between frames approaches the particle 
spacing. 

/From the particle tracks we determine average parti- 
cle velocities V{y) and RMS velocity fluctuations perpen- 
dicular to the flow direction (8V x (y)) and parallel to it 
(6Vy(y)) as a function of distance y from the rotating in- 
ner cylinder. The position resolution of < 0.1 pixel yields 
a resolution of particle velocities and fluctuations of bet- 
ter than 0.1 pixel/frame. The upper limit for measurable 
velocities is given by the tracking routine, which requires 
that the maximum displacement be smaller than one par- 
ticle diameter. Larger particle images yield a larger ve- 
locity range, but the velocity profile may not improve 
since fewer particles can be tracked in a single image. 
A mean particle size of about 20 pixels gives sufficient 
dynamic range for the velocities with good statistics. 

IV. EXPERIMENTAL RESULTS 

A. Particle Dynamics 

The behavior of the inner cylinder is found to be very 
similar to the dynamics of a rough plate sliding across a 
granular layer |2{J . At low shear rates the motion of the 
inner cylinder is intermittent with short, rapid slips, and 
long periods of sticking. At sufficiently high shear rates 
or with a stiff connection between motor and cylinder, 
steady motion of the inner cylinder is observed. Air flow 
reduces the shear forces, as already noted by Tardos et 
al. ||, and it also suppresses stick-slip motion. 
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FIG. 2. Mean particle velocity (normalized by the shear 
rate) as a function of distance from the inner cylinder (in 
particle diameters) for glass spheres. The respective shear 
rates U (in Hz) are : 0.004 (hexagons), 0.04 (squares), 0.01 
(open triangles), 0.4 (crosses). The solid triangles show the 
velocity profile at U = 0.01 Hz without air flow. The normal- 
ized mean velocity profile is independent of shear rate and 
shear dynamics (intermittent or steady motion). The dashed 
line is the solution of (|2^), with S = 4.7d, y w = 2.8d, and 
a = 0.4 ( see text for details). 

The velocity profile V(y), when normalized by the 
shear rate, is roughly independent of shear rate, as shown 
in Fig. ^. The mean velocity profile without air flow in 
the stick-slip regime (solid triangles) is essentially the 
same as that for steady shearing. The dashed line shows 
the theoretical prediction from our hydrodynamic model, 
which will be discussed in section [V| 

Fig. U shows the perpendicular RMS velocity fluctua- 
tions SV y , which have not been previously measured in a 
3D system to our knowledge. When data taken at differ- 
ent shear rates U are normalized to the same magnitude 
at a distance of 3 particle diameters away from the wall, 
the fluctuations follow the same profile, independent of 
shear rate and independent of the presence or absence 
of stick-slip motion. The velocity fluctuations decrease 
roughly exponentially far from the inner cylinder but fall 
off more slowly with y than does the average velocity. 
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FIG. 3. RMS velocity fluctuations perpendicular to the 
shear direction. Fluctuations decrease roughly exponentially 
far from the inner cylinder, but more slowly than the mean 
flow. The rms fluctuations are rescaled (shifted vertically) 
such that all experimental points are forced to agree at y — 3d, 
where d is the bead diameter. The dashed line is the theo- 
retical result (see text), with a decay length 8 — 4.7rf and a 
boundary position y w = 2.8d. Measurements are made on 
glass spheres, with symbols as in Fig. 2. 

We determine fluctuations by averaging the velocity of 
a small section of the image over a long time and then cal- 
culating deviations of individual particle velocities from 
the mean of that section. The measured parallel fluctua- 
tion amplitude therefore includes fluctuations of the flow 
speed, which can be caused by the soft spring connec- 
tion to the motor. In order to compare velocity fluctua- 
tions during stick-slip motion and steady shearing and to 
compare the measurements to the hydrodynamic model 
of a steady state flow, we show only the perpendicular 
fluctuations. We note that even during steady shearing, 
parallel fluctuations are larger than perpendicular fluc- 
tuations, but their ratio remains roughly independent of 
y, as shown in Fig. 5 below. 

In principle the density profile could also be measured 
using the tracking algorithm, by counting the average 
number of tracked particles as a function of position y. 
In pratice however no quantitative results could be ob- 
tained because of the limitation of the tracking method to 
resolve particle positions in the third dimension (e.g., for 
low densities, particles from lower layers are also counted: 
the density is thus quantitatively overestimated, while 
this does not affect the mean and RMS velocity profiles). 
Qualitatively, the measured density profile increases with 
the radial coordinate toward a limiting value at large 
distances (not shown). The density close to the moving 
boundary is measured to be up to 40% below its lim- 
iting value, depending on shear rate and airflow. This 



is in agreement with other numerical and experimental 
observations 34 1. 

We now examine the particle behavior near the bound- 
aries in more detail. Close to the inner cylinder it is 
difficult to distinguish wall particles from particles that 
move close to the wall. Because we image the surface 
from a slight angle and because the height of particles 
fluctuates slightly, the boundary between wall particles 
and sheared particles fluctuates. We have examined the 
boundary conditions with both steel spheres and rough 
ceramic spheres, which allow us to distinguish the par- 
ticles from the layer of rough glass beads glued to the 
inner cylinder. We find that the granular temperature 
has an approximately constant value in a region about 3 
particle diameters wide near the inner wall (see Fig. 3). 

Since particles barely move close to the stationary 
outer cylinder, the granular temperature at large y is 
examined at a lower frame rate than is necessary at po- 
sitions close to the rapidly moving inner cylinder. This 
yields better statistics for fluctuations in particle motion. 
Fig. |J shows the mean velocity and velocity fluctuations 
close to the outer wall. The velocity decreases roughly 
exponentially up to roughly two particle diameters from 
the outer wall. The velocity fluctuations also decrease 
approximately exponentially in this section. 




FIG. 4. Velocity(solid circles) and fluctuation profile (open 
circles) close to the stationary outer cylinder. Each profile is 
normalized by the shear rate. 

As discussed in the previous section, an upward airflow 
reduces the effective pressure, while a downward airflow 
increases the effective pressure within the material. This 
allows us to measure the pressure dependence of the ve- 
locity and fluctuation profile. Fig. |s| shows three exper- 
iments at different air flow rates (i.e. pressure). In or- 
der to avoid stick-slip motion without airflow and with 
downward airflow, the motor is connected rigidly to the 
cylinder for these experiments. The crossing points of 
the temperature and velocity profiles in Figs. 4,5 are 
different because the mean velocity was different. 
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Neither the velocity profile nor the profile of RMS fluc- 
tuations (the granular temperature) change with pressure 
over the range of pressures accessible with this method. 
This is consistent with our hydrodynamic model as de- 
scribed in section |y|. Note that the RMS fluctuations 
parallel to the shear direction are larger by a factor of 
roughly 1.3, even though the mean velocity cannot fluctu- 
ate due to the rigid connection between motor and cylin- 
der. This anisotropy has been observed previously fl30| . 
It may be connected to an anisotropy in pressure in a 
sheared granular system plp^ ]. 
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FIG. 5. Velocity profile V{y) (solid symbols) and RMS fluc- 
tuation profile SV(y) (open symbols) both perpendicular (no 
lines) and parallel (lines) to the shear direction at different 
effective pressures controlled by air flow. Measurements on 
glass spheres with upward airflow (triangles pointing up), 
no airflow (circles), and downward airflow (triangles point- 
ing down). The profiles do not depend significantly on the 
pressure. 

The velocity and fluctuation profiles are independent 
of shear rate and shear dynamics, but can vary with the 
material that is sheared. Our model contains three pa- 
rameters that determine the velocity and fluctuation pro- 
files, as discussed in section These parameters may 
depend on material properties in a non-trivial way. In 
Fig. H we compare three very different materials (steel 
spheres, glass spheres, rough porous ceramic particles). 
For all materials, the mean velocity decays roughly expo- 
nentially far from the shear boundary. The characteristic 
length of that decay is between 1.5 and 2 particle diame- 
ters. The velocity of glass spheres decreases more slowly 
away from the sheared cylinder than does that of rough 
ceramic particles. Steel spheres are about twice as large 
as the glass spheres coating the surface. This leads to 
significant slip at the boundary. The RMS fluctuations 
decay significantly more slowly than the mean velocity 
for all materials (not shown). 




Particle diameters 

FIG. 6. Comparison of V(y) for smooth glass spheres (cir- 
cles), rough porous ceramic particles (triangles), and steel 
sphere mixtures (squares) at U = 0.04Hz. 

We have also done experiments in a narrow gap geome- 
try, where the shear region is only 4 — 5 particle diameters 
thick. In this case the velocity profile is linear and the 
temperature is roughly constant across the cell, as shown 
in figure [?]. 
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FIG. 7. Velocity profiles (Solid lines) and RMS fluctua- 
tion profiles (dashed lines) in a narrow gap geometry with 
gap width of 4 — 5 particle diameters. The velocity profile is 
linear and the RMS fluctuation profile is constant. Measure- 
ments are made on glass spheres with U = 0.1Hz (circles) and 
U = 0.02Hz (triangles). 

Fluctuations parallel to the shear are somewhat larger 
than perpendicular fluctuations. The velocity profile is 
again independent of shear rate. The small gap result 
is consistent with the hydrodynamic model of section [v|: 
The uniform RMS fluctuations across the gap are accom- 
panied by a linear velocity gradient. 
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B. Shear forces 

The shear stress is found to be roughly independent of 
shear rate but to decrease roughly linearly with increas- 
ing upward air flow as shown in figure ||. The dependence 
on air flow is consistent with the results by Tardos et 
al. Q described in section H. Previous experiments |3^ 1 
showed that the shear stress is directly proportional to 
the pressure inside the granular material. We can there- 
fore assume that air flow decreases the pressure roughly 
linearly with increasing flow rate. 

When there is no air flow, we also find that the shear 
stress is roughly independent of shear rate even though 
stick-slip motion is observed. This indicates that some 
velocity weakening (i.e., a decrease in shear force with 
increasing velocity) must occur. The mean shear stress 
with air flow is a factor of four smaller than the mean 
shear stress without air flow. For a small gap, the 
shear force increases with decreasing shear rate, even- 
tually leading to jamming of the inner cylinder below a 
threshold shear rate. 
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FIG. 8. Mean shear stress vs air flow rate at shear rates 
of 0.004 Hz, 0.01 Hz, 0.04 Hz, and 0.1 Hz. The shear stress 
decreases approximately linearly with increasing air flow. The 
shear stress is independent of shear rate at most air flows. 
(Glass spheres) 

The air flow at which the transition from stick-slip mo- 
tion to steady shearing is observed is shown as a function 
of shear rate in Fig. [|. We determine the transition from 
stick-slip motion to steady sliding motion from the emer- 
gence of a peak at v = in the probability distribution of 
the shear velocity. The critical air flow decreases roughly 
linearly with shear rate. 
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FIG. 9. Air flow at which the transition from stick-slip mo- 
tion to steady shearing is observed as a function of the shear 
rate. The critical air flow decreases roughly linearly with 
shear rate. (Glass spheres) 



V. THEORY 

A. Hydrodynamic Model and main assumptions 

Here, we propose a hydrodynamic model for granu- 
lar flow to describe the data presented in the previous 
section. By hydrodynamic, we mean that a local, mean 
relationship is assumed to hold between shear stress and 
shear rate, in contrast to recent approaches which ad- 
vocate nonlocal relationships (see e.g., |2C|1 , and yj and 
references therein). The grains are assumed to behave 
like inelastic hard spheres, with a diameter d and coef- 
ficient of restitution e. Wc therefore neglect any friction 
force between grains. Within this model, collisions be- 
tween grains are instantaneous. The inelasticity coeffi- 
cient e is moreover assumed to be independent of the 
relative velocity of the two colliding particles. Our phi- 
losophy is different from more sophisticated approaches 
(such as the simulations of Ref. |34j), where the micro- 
scopic model is taken to be as realistic as possible. Here 
we deliberately choose the simplest model: our goal is to 
show that despite its simplicity, a hydrodynamic model 
leads to many non-trivial results which are usually at- 
tributed in the literature to more sophisticated ingredi- 
ents of granular flow. We check its validity by comparing 
a posteriori our findings to the experimental results. As 
we demonstrate below, most of the experimental prop- 
erties of sheared granular flow can be explained by this 
hydrodynamic, locally Newtonian description of the ma- 
terial. 

We consider a simplified Couette geometry: the gran- 
ular material is confined between two parallel walls, sep- 
arated by a distance H. The x axis is along the walls, 
while the y axis is perpendicular to it. The bottom wall, 
placed at x — 0, is assumed to move at a velocity U along 
the x direction and the top wall stays at rest. 
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We start with the equations of hydrodynamics for the 
inelastic hard sphere fluid. Following common practice, 
we identify the granular temperature with m[(Sy] 2 , where 
SV is one component (x or y) of the RMS velocity fluc- 
tuation previously defined in section [II, and m is the 
particle mass. There is a small ambiguity in this defi- 
nition for a sheared granular material since the x and y 
components of [5V] 2 (i.e. the components parallel and 
perpendicular to the flow) differ slightly. However, as 
may be observed in Fig. |5|, they differ only by a scaling 
factor close to unity. We shall therefore disregard this 
anisotropy in our theoretical treatment. For the inelastic 
hard sphere fluid of local particle density p, the equa- 
tions of hydrodynamics for the mean velocity field V and 
temperature T can be written as || 



'dV Tr „ T \ „ 
mp\ — + V ■ W I = — V • o~ 



mp 



at 

fdT 



(la) 



+ V • VTJ = -V • Q - cr : k - eT, (lb) 



where : means contraction of the two tensors. In these 
equations, k is the symmetrized velocity-gradient tensor 



K>a,0 = g {d<*Vp + dpV a ) , 



(2) 



a is the pressure tensor, Q is the heat flux, and e is the 
temperature loss rate per unit volume. As in a Newto- 
nian fluid, we assume a linear, local relationship between 
fluxes and forces. We thus write the pressure tensor as 



a = PI - 2r) (« - V • VI) , 



(3) 



where P is the pressure, r\ is the shear viscosity and I the 
unit tensor. In a similar way, we assume Fourier's law for 
the heat flux, 



Q 



-AVT, 



(4) 



with A the thermal conductivity. These equations are 
completed by the equation of state of the material, in 
the form P = P(p,T). 

In the simplified planar shear cell geometry, the mean 
flow is a function of y only and parallel to the x direc- 
tion: V = V(y)e x with e x the unit vector in the x direc- 
tion. ^From the momentum transport equation (|la| ) in 
the steady state one then expects 



d 



(5) 
(6) 



This shows that both the shear stress a xy (=o~ yx ) and 
the pressure P = a yy are independent of y. 



B. Equation of state and high-density expressions 
for the transport coefficients 

In order to solve the hydrodynamic equations (]l|), ex- 
plicit expressions for the transport coefficients and the 
equation of state in terms of density and temperature 
are needed. 

When the system is at rest, the density p in the granu- 
lar material is roughly given by the random close packing 
(RCP) density p c (the granular system does not to crys- 
tallize because of a slight polydispersity of the beads). 
Experimentally shear rate is found to dilate the sys- 
tem ||. More precisely the density decrease is larger 
close to the boundaries, where the shear rate is larger, 
than far from the moving wall, where the velocity goes 
to zero [§jB4j]- As already quoted in Sec. III., our ex- 
perimental results for the density do qualitatively agree 
with this observation. Consequently a broad range of 
densities is explored in the sheared system, going from a 
value slightly below the RCP density far from the mov- 
ing wall, down to a density up to 40% less than RCP 
near the moving wall at high shear rates. This makes 
the problem much more difficult than in standard fluids, 
where the density remains constant over the cell. Here 
the functional dependence of the transport coefficients on 
density is an important ingredient of the theory since it 
does affect flow properties. 



1. Equation of state 

For an inelastic hard sphere system, the equation of 
state can be written in terms of g(d) the pair correlation 
function at contact (d being the diameter) |35],[56) , in the 
form 



P = pT 



l + (l + e)-pd 3 g(d) 



(7) 



(p is the local numerical density, i.e. the number of par- 
ticles per unit volume). For both dilute and moderately 
dense systems (pd 3 ~ 1), g(d) is accurately described by 
the Carnahan-Starling formula [[35|. However since the 
density of the granular material in the shear cell is close 
to random close packing (RCP) where g(d) diverges, an 
alternative expression for g(d) is usually assumed J37|: 



9(d) 



i - pIpc 



(8) 



with p c the density at RCP. The equation of state, Eqs. 
([?])- (||), then takes the following approximate form in the 
high density limit 



P = Po 



(1 - ■£-) 

^ Pa I 



T, 



(9) 



with p = (1 + e)^p^d 3 . The equation of state thus gives 
a natural link between temperature and density: since 
the pressure is independent of y, one obtains 







1 



£_ 

Pc 



PO 
P 



T. 



(10) 



This equation shows that the regions of small temper- 
ature correspond to high density regions, p ~ p c , and 
vice- versa. 



2. Transport coefficients 

Expressions for the transport coefficients of the inelas- 
tic hard-sphere model have been computed from the En- 
skog equation |^,36||. This kinetic equation takes the ef- 
fects of excluded volume into account, but it neglects 
any correlation between the velocities of colliding parti- 
cles. This approach yields expressions for the transport 
coefficients in terms of the density, temperature, and pair 
correlation function at contact. We refer to Refs. [p|j36| 
for explicit expressions. We note that in the high den- 
sity limit, g(d) becomes very large and the dependence 
of the Enskog transport coefficients on density mainly 
come from the terms proportional to g(d). In this case, 
the transport coefficients reduce to the generic forms 



r) E (p, T)~r] - 



,1/2 



<P 



fJ 



(d) T 1/2 



\ E (p,T) ~ Ao^r- g(d) T 1 ' 2 



m 1 ' 2 d 2 
1 

m l l 2 d 



e E (p,T) ~ e p l j2 g{d) T 1/2 , 



(11) 



where m is the mass of the particles and 770, Ao, eo are 
dimensionless coefficients that depend only slightly on 
density in the high-density limit. If p ~ p c as discussed 
above, these coefficients can be taken as numerical con- 
stants. Finally, we mention that eo is proportional to 
1 — e 2 , where e is the coefficient of restitution. Thus 
eo = in the purely elastic case, as expected. 

The full expressions for the transport coefficients ob- 
tained within the Enskog kinetic theory have been found 
in simulations to be correct for small and intermediate 
densities [p4| . However, as mentioned above, the den- 
sity of the flowing material in the shear cell is higher 
and close to the RCP density. In this limit, some of the 
Enskog expressions for the transport coefficients may no 
longer be valid, mainly because correlations between col- 
liding particles and collective phenomena, which are not 
included in the Enskog theory, then play an important 
role. The reason is quite intuitive: at high density close 
to random close packing, a particle can move substan- 
tially (over a distance of the order of its diameter) only 
if its neighbors move coherently. Only collective motion 
is therefore possible. It has been found in Molecular Dy- 
namics simulations of the hard sphere model [^4j that 
these correlations only affect the shear viscosity and self 
diffusion coefficient, which depart from their Enskog ap- 
proximation at high densities. On the other hand, the 
thermal conductivity has been found to be well described 
by the Enskog expression up to very high densities. Such 



behavior is in fact expected since transport of energy 
does not require motion of particles over large distances 
(only "rattling" around the mean position is involved in 
energy transport). Similar conclusions were reached by 
Leutheusser |2j| on the basis of a mode-coupling calcu- 
lation for the clastic hard-sphere model. 

As a result of these considerations, the thermal conduc- 
tivity A and loss rate e are expected to keep their Enskog 
expressions \ E and ce, as given in Eqs. (pd|), over the 
whole range of densities (between intermediate densities 
up to RCP): 

A(p,T)~A - ' ' 



e(p,T)~ e - 



1 



1 



JTl/2 



lV2d 



1 



T 



1/2 



(12) 



On the other hand, a crossover is expected for the 
shear viscosity between its Enskog approximation, in 
Eq. ([ll]), for intermediate densities, towards an asymp- 
totic stronger divergence as a function of density very 
close to RCP. Such a crossover is indeed found in Molec- 
ular Dynamics calculations of the self diffusion coefficient 
in a monodisperse hard sphere system close to freez- 
ing |58j . By analogy with the behavior of supercooled liq- 
uids above the glass transition |^| , we shall assume that 
very close to RCP, the viscosity diverges algebraically as 
a function of density p near p c : 

^D-^L^i-v*. (i3) 



('-*) 



Here 771 is a dimensionless numerical constant. At this 
point the exponent is a phenomenological parameter 
in the theory. Since the viscosity is expected to diverge 
more strongly than its Enskog expression, we expect [3 to 
be larger than one. We shall discuss in the next section 
how the form proposed in Eq. ( |ilj| ) compares with the 
experimental results. 

We emphasize again that such an algebraic divergence 
of the viscosity is expected in supercooled liquids close 
to the glass transition [j59). More precisely, an algebraic 
divergence is associated with the existence of cooperative 
interactions which predominate over thermall y a ctivated 
processes |p9fl . The functional form in Eq. ( |l3|) is pre- 
dicted by mode-coupling calculations of the dynamics of 
supercooled liquids. In the case of 3D hard sphere, the 
latter approach yields an exponent (3 = 2.58 pQ |. 

To our knowledge, there is no expression for the viscos- 
ity over a full range of densities that would make explicit 
the cross-over between the Enskog expression at inter- 
mediate densities and the asymptotic expression at very 
high densities. In order to avoid restrictive assumptions 
about the general form of the density dependence of the 
viscosity, we shall write in general 



V(p,T) = fj(p)^T^, 



(14) 
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where the function rj(p) has the two limiting forms: 

no 



p~d A < p c r)(p) 



P ~ Pc 



v(p) 



m 



(15) 



As we shall show in the following sections, the knowledge 
of these two limiting behaviors is sufficient to obtain a 
quantitative description of the flow and a qualitative pic- 
ture for the shear forces. 

Combining Eqs. @ and (jl^) allows us to write the 
transport coefficients in terms of the pressure and tem- 
perature as 



n(p,Ty 



,V2 



d 2 \P J 

\(p,T)~\ ^ P 
m 1 ' 



(PqT\ T l/2 



e{p,T)~e — 

m - 1 



m 1 / 2 d 2 poT 1 / 2 
1 P 



m 



l 2 d T 1 / 2 ' 



(16) 
(17) 
(18) 



(Note that in order to improve readability, we dropped a 
numerical constant prefactor p c / Po m the expression for 
e, which amounts to a rescaling of the numerical prefactor 



C. RMS and Mean flow velocity profiles 

In this section we compute the granular temperature 
profile, T(y) and mean- velocity profile, V(y). In the sta- 
tionary Couette geometry, the hydrodynamic equation 
for T(y) is found to reduce to 



— \(p T)—T 

dy ' dy 



Oxyi - t{p,T)T = 0, 



and 



o X y = nip, T )i = const, 



(19) 



(20) 



where 7 = dv ^ v ^ 

dy 



is the shear rate. Using the expres- 



sion for the shear stress, a xy = r)(p, T)j, the second 



Jn(p, T )- 



xyl — "xyl 

dvy/dx is negative 



term of Eq. ( |i9| ) can be rewritten as a 
Note that a xy is negative because 7 
throughout the sample. 

Both equations involve the explicit form of the density 
dependence of the viscosity, ?y(p, T), for which no explicit 
functional form is available to us. At first sight, it would 
seem hopeless to obtain a full expression of the temper- 
ature and velocity profiles. This is not the case. As 
we shall show, a simple phenomenological picture, which 
emerges from the asymptotic forms of r)(p, T) as given in 
Eq. (Qq) , allows one to overcome this problem and to ob- 
tain tractable expressions for the temperature and mean 
velocity. 



We proceed in two steps. First, we describe this phe- 
nomenological "two region" picture and obtain expres- 
sions for the velocity and temperature profile. Then in 
a second step we come back to a more general but for- 
mal solution of Eqs. ( |l9| ) and (^0|) (in the next section). 
This general discussion allows us to discuss the velocity 
dependences of the shear forces. 

We start with the discussion of the temperature pro- 
file. The role of the the nonlinear term a xy j in Eq. © 
is in fact physically quite simple to understand. It merely 
acts as a source ("heating") term for the fluctuations: it 
is through this nonlinear contribution that the flow cre- 
ates the fluctuations which in fine couple back to the 
mean flow. However, this source term is only effective 
close to the wall as we show below. This simplifies con- 
siderably the picture for the creation and transport of 
temperature. Basically, two regions can be treated sep- 
arately: close to the moving boundary, fluctuations are 
"created" through the nonlinear coupling to the flow; far 
from the boundaries, inelastic dissipation dominates over 
nonlinear heating, and the local temperature profile is de- 
termined entirely by the diffusion and heat loss terms of 
Eq. ©. 

This separation can be understood on the basis of the 
asymptotic behavior of the viscosity discussed above. In- 
deed, it is easy to show that far from the wall the non- 
linear term a xy j goes to zero faster than the dissipa- 
tion term e(p,T)T in Eq. @. Away fr om the moving 
boundary the temperature goes to zero and the density 
goes accordingly towards RCP. From Eqs. ([l5]) and J18] ) 
it follows that the nonlinear term behaves in this region 
like <j 2 y /r](p,T) oc T^ -1 )/ 2 , while the dissipative term 

in this region scale with temperature as e(p, T)T oc T 1 ! 2 
(where Eq. ( |l8l ) has been used). Since we anticipated 
that the exponent j3 is larger than one (around 1.75 as 
found experimentally, see below), the nonlinear term de- 
cays more strongly towards zero than the dissipative one. 
On the other hand, the nonlinear term is relevant close 
to the moving boundary where the shear rate is large (or 
equivalently the density is smaller). We note that if the 
exponent had been equal to 1, both previous terms 
would have been comparable and then the picture of two 
separate regions would not have been appropriate. 

Close to the walls, the full nonlinear equation includ- 
ing the "heating" term should be solved. This is in fact 
not necessary to make simple predictions for the RMS 
and mean velocity profiles. An inspection of the experi- 
mental profiles, as presented in Fig. 0, shows that in the 
region close to the moving boundary, the temperature is 
roughly constant over a small layer, several particle di- 
ameters thick, with a thickness which can be considered 
in a first approximation to be independent of the shear- 
ing velocity U. As a first step, we therefore assume in 
an ad hoc way that for distances y smaller than a cut-off 
distance y w , the temperature T is constant, T(y) = Jo- 
in this pictorial view, the boundary layer corresponds to 
the region where the nonlinear term is important. At 



ff 



this stage, the parameter To merely plays the role of a 
normalizing constant. We shall come back to this point 
in much more detail in the discussion of the shear forces 
because, while the precise value of To does not influence 
the velocity profile, it does strongly affect the prediction 
for the shear forces. 

We are now in a position to obtain an explicit expres- 
sion for the temperature profile. For distances larger than 
y w , the transport equation for the temperature, Eq. (19), 
reduces to 



-Uo — 

rill 1 ml/' 



P d 



dy \ "mMH 2 poT 1 / 2 dy 



T-e - 



1 



iV 2 d 



p T l/2 = Q) 

(21) 



where the high-density expressions for the transport co- 
efficients, Eqs. (|lg|), has been used. 

As shown previously, the pressure P is independent of 
y, so one may rewrite Eq. (|2l| ) as 

_g_ rl/2 _ l_ T ^ = Q) (22) 

where 5 has the dimension of a length and is defined as 
c2 2Aq 



eoPo d 



(23) 



The parameters Ao, and eo are just numerical constants 
in the high-density regime of interest to us. Therefore, 
one expects S to be of the order of a few particles diam- 
eters. Note that since eo is proportional to 1 — e 2 (with 
e the restitution coefficient), the decay length 5 goes to 
infinity when the system become perfectly elastic, as one 
would expect. This equation has to be completed by 
boundary conditions for the temperature at both walls. 
At the moving wall, we set T(y) = Tq for y < y w , as 
discussed above. At the wall at rest, detailed experi- 
ments, as shown in Fig. ^, show that the temperature 
profile is compatible with a vanishing heat flux condition 
dT/dy = 0. (In general, one expects the boundary con- 
dition for the temperature to relate the heat flux at the 
boundary to the product of the interface (Kapitza) resis- 
tance and the temperature jump: Here we just assume 
that the Kapitza resistance is very small. We shall come 
back to this point in the next section.) 

The solution of Eq. (^) with these boundary condi- 
tions is 



V < Vw T(y) = T 
y w <y T^ 2 (y) = T ( ) 



coshf^ 

1/2 V S 



cosh 



H-y u 



(24) 
(25) 



In this equation, H is the thickness of the shear cell. As 
shown in Fig. ^ this result is in good agreement with the 
experimental data. We find that S = 4.7c? and y w = 2.8d 
for glass spheres, but these parameters depend somewhat 
on the material properties. 



The mean velocity profile can be obtained along the 
same lines from the temperature profile of Eq. (|25|). 
Eq. (20) gives unambiguously the shear rate 7 in terms 



of the temperature profile T(y), since the y dependence 
of the density is already contained in the temperature 
through the equation of state, Eq. ([io|). It seems at first 
sight difficult to make an explicit prediction for the veloc- 
ity profile, since the explicit expression of the viscosity is 
not known over the full range of densities. However, the 
"two regions" picture, which emerged in the discussion of 
the temperature profile, is relevant for the mean- velocity 
profile as well. Far from the wall (i.e. for y > y w ), 
the temperature decays to zero and the density goes ac- 
cordingly to the RCP limit (as shown in Eq. (10)). In 
this region, one thus expects an "anomalous" scaling for 
the density dependence of the viscosity, according to Eq. 
(|l5|). Thus one finds for y > y w 



,1/2 



Vo- 



d 2 



, — y 



( 



7 



(26) 



It is interesting to note that this relation yields a scaling 
relationship between the temperature T and 7 that can 
be tested experimentally. The previous equation can be 
rewritten in the form 



(0 



J 1 ! 



(27) 



This power law relationship has been checked experimen- 
tally, as shown in Fig. [n], where the local RMS velocity 
profile SV = [T/m) 1 / 2 is plotted versus the local velocity 
gradient dV x /dy on logarithmic scales. 
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FIG. 10. Connection between the local RMS velocity fluc- 
tuations and local shear rate (same symbols as for Fig. |2|). 
Local fluctuations are found to increase approximately as a 
power law of the local velocity gradient, with a power of 0.4 
(dashed line). 
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tern. We have a 
Howell et al. fET 



In this plot a scaling exponent a = 0.4 is obtained for 
the glass bead system. This allows for an experimental 
determination of the exponent (3 for the "anomalous" di- 
vergence of the viscosity close to RCP. According to Eq. 
(p7|), a is identified in the theory as a — 1/(2/3 — 1), 
so that (3 = 1.75. The fact that this exponent clearly 
exceeds unity confirms the collective character of the dy- 
namics in the granular material. This value of the ex- 
ponent is slightly smaller than the previously mentioned 
result for (3 = 2.58 obtained within the mode coupling 
theory applied to the case of hard spheres pOfl . A pos- 
sible reason for this difference could lie in the rotational 
degrees of freedom of the beads in the experimental sys- 
tem which are absent from the mode coupling estimate. 
By increasing the number of degrees of freedom an ef- 
fective decrease of the shear viscosity could be obtained. 
We leave this question for further investigation. 

We emphasize that the observation of the scaling rela- 
tionship between 7 and T 1 / 2 is not restricted to our sys- 
iplied the same procedure to the data of 
which were obtained in a 2D system of 
photoelastic disks : for the 14 different densities studied 
in this work, a similar scaling relationship between the 
local 5Vrms and 7 is measured. An exponent a ~ 0.5 is 
obtained, yielding f3 ~ 1.5 for that experiment. Moreover 
both the mean and fluctuation profiles obtained by How- 
ell et al. can be fitted within the present model. This 
observation seems to indicate that the coupling between 
mean and RMS fluctuations which is the essence of our 
model is a general feature of the underlying dynamics of 
the granular material. 

Once the exponent p is determined, the velocity profile 
can be obtained by integration. In the boundary layer for 
the temperature (y < y w ), the temperature is assumed 
to be a constant; Eq. ( p0| ) then shows that the shear rate 
is constant, yielding a linear velocity profile. This ob- 
servation is in agreement with the experimental velocity 
profiles, as can be seen by a careful inspection of Fig. g 
for y < y w . On the other hand for y > y w the velocity 
profile is obtained by integrating Eq. ( p7| ) together with 
the solution j25|) for the temperature profile. The full ve- 
locity profile is obtained by matching the two solutions 
at y w (where we assume continuity of velocity and shear 
rate). The solution obeying no-slip boundary conditions 
at both walls takes the form 



V x (y) = U\l 



lody' M) 



Io H dy' <t>{y>) 
where the function <f>(y) is defined as 

y < y w <t>(y) = 1 

cosh(^) " 



y>yv 



Hy) 



cosh(^) 



2/3-1 



(28) 

(29) 
(30) 



When the value (3 = 1.75 previously obtained for the 
exponent (3 is used, the theoretical velocity profile defined 



in Eqs. ( |28|J30| ) is easily integrated numerically. The re- 
sult is plotted in Fig. |2| together with the experimental 
results. As can be seen, good agreement is obtained with 
the theory. 

Although the two region picture seems to be quite suc- 
cessful in describing the mean and RMS velocity profiles, 
it would be more satisfactory to have a full expression 
of the density dependence of the viscosity over the full 
range of densities in order to integrate Eqs. ( |l^ ) and ( |20| ) 
explicitly. The two region picture is to be considered 
merely as a simple and physically sound way of dealing 
with our ignorance of the explicit functional dependence 
of the viscosity on density. Its success indicates that the 
velocity profile does not depend crucially on the details 
of this relationship. 



D. Towards "Solid like" Shear Forces 



1. General discussion 

In his pioneering work, Bagnold measured the shear 
force in a Couette cell and found a quadratic increase 
of the shear force as a function of the shearing velocity, 
a xy tx U 2 Q). Using kinetic arguments he proposed a 
phenomenological relationship between shear stress and 
shear rate in the form a xy oc 7 2 , which can be under- 
stood by assuming that the collision frequency between 
granular particles is fixed by the shear rate 7 itself. Sub- 
sequent experiments were performed @] which showed 
that the functional dependence of the shear force as a 
function of the shearing velocity does depend on whether 
the system is allowed to dilate or not: when the system 
is sheared at constant volume, the shear force is found to 
be proportional to U 2 (i.e. the Bagnold scaling); when 
the system is allowed to dilate, the shear force is found to 
be independent of U. In the present geometry, the mate- 
rial is allowed to dilate and a solid-like friction force,i.e. 
approximate independence of U, is found. 

It is tempting to ascribe the solid like behavior to a 
solid friction force existing between granular particles 
at contact. On the other hand it is difficult to recon- 
cile this behavior with constant-volume Bagnold scaling, 
a xy oc U 2 , whose origin is purely kinetic. Here we shall 
show that both behaviors can be accounted for within 
the hydrodynamic model we have introduced. 

• Constant volume results: 

It is evident that the Bagnold scaling, a xy oc U 2 , 
originates from a general dimensional argument, 
whatever the local relationship between shear stress 
and shear rate is: in other words, the global scal- 
ing a xy oc U 2 obtained in a constant volume ex- 
periment does not necessary imply that the local 
relationship a xy oc 7 2 . The reason is very simple. 
Let's consider a system of N inelastic hard spheres 
contained in a plane Couette cell of (fixed) volume 
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V. The bottom boundary is moving at a constant 
velocity U, while the upper boundary stays at rest. 

In this case, the only microscopic velocity scale en- 
tering the dynamics is U, the velocity of the mov- 
ing wall. In other words, all dynamical quantities 
can be written in terms of the shear velocity U 
only, which fixes the only time scale in the problem. 
^From the analysis of the full TV-body dynamics, 
one thus concludes in general that 



Vx{y) 
u 

T(y) 

mil 2 
P 

O xy 



= f 



y N H 
d' V' ~d 
y N H 
d'V'~d 



N 



N H 



N 



N H 



v a ° {V ~d 



mU z 



raU 1 



(31) 
(32) 



where /, g, p, o~q are dimensionless functions; H is 
the cell width and d the bead diameter. 

In this context, the Bagnold scaling results natu- 
rally from the last equality in (|32|), whatever the 
relationship between o~ xy and 7 is. The only impor- 
tant condition is that the microscopic dynamics of 
the granular material only involve binary collisions. 

• Constant pressure results: 

Experiments performed in a constant pressure sit- 
uation are more subtle to handle. The introduction 
of a constant pressure in the system introduces a 
new time scale, namely (m/ Pd) 1 ^ 2 , which will com- 
pete with the time scale associated with the shear- 
ing velocity, d/U. Thus different regimes might be 
found depending on whether the time scale associ- 
ated with the shear rate is larger or smaller than 
the one associated with pressure. 

We show in the following that in the situation of con- 
stant pressure, a shear force independent of shear rate 
is obtained at large shearing velocity, while a velocity 
weakening regime is found at small velocities. 



2. Shear forces at constant pressure from the 
model: qualitative 



In the previous section the experimental mean and 
RMS velocity profiles were reproduced by using a sim- 
ple " two regions" picture, which allowed us to overcome 
our lack of knowledge of the full functional dependence 
of the viscosity on the density. However the shear force 
originates in the interaction between the wall and the 
fluid particles adjacent to it. Therefore one needs a good 
description of the region close to the wall in order to 



get valuable information about the shear forces. This re- 
quires a careful analysis of the nonlinear "heating" term 
in the equation for the temperature, Eq. (|l9|), which we 
perform in this section. 

We first present an approximate, qualitative argu- 
ment that helps explain physically what quantities con- 
trol the friction force in the shearing experiment. We 
then present the full approach, which clarifies the crude 
approximations made in the qualitative arguments. 

The shear force is defined in the hydrodynamics model 
as the viscosity times the shear rate, both evaluated at the 
boundary. Using the expression for the viscosity, Eq. jl4|), 
this yields: 



V 



f PqTq 



,1/2 



-T, 



1/2. 



7 lo- 



We need a second equation giving the temperature at the 
boundary T . In the "two regions" picture we introduced 
above, the temperature is created near the boundary and 
then transported to the rest of the system In this picture, 
the temperature T merely comes from a balance between 
the dissipation rate close to the moving boundary — (Tq 
and the "heating term", a xy j\o. Together this gives 



O-xy"/\0 



£ p T l/2 



(34) 



where we have used the expression in Eq. ([18]) for the 
dissipation rate e. 

The shear rate 7I0 can be eliminated from the two pre- 
vious equations, Eqs. (J33|) and (34), to give 



^ xy 
~P~ 



{pod' 



PoT \ p T 



P 



P 
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(35) 



Now, depending on the ratio pqTq/P, two limiting be- 
haviors are obtained 

(i) If PqTq/P is large, then the viscosity takes its in- 
termediate density expression in Eq. (|l5|), fjippTo/P) ^ 
ijqP/ PqTq. In this limit, one then obtains from Eq. ( p5|) 



'xy 



PqP 



(36) 



with pq — eQ-rjo/pod 3 a dimensionless constant. In other 
words, the shear force is expected to be independent of 
the shear rate and proportional to pressure, as is usually 
found in solid friction. 

(ii) If PqTq jP is small, then the viscosity takes its high 
density expression in Eq. Jl5|), ^(pqTq/P) ~ r)i(P / PqTq) 13 
and Eq. (133) gives 



^ xy 
~P~ 



( P 



p d 3 \pqTq 



0-V 



1/2 



(37) 



The temperature Tq enters explicitly in this formula. It 
has to be obtained as a function of the shearing veloc- 
ity U. This can be accomplished using the energy bal- 
ance discussed above, which results in Eq. (^3). In this 
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equation, one expects the shear rate at the boundary to 
be roughly given by -y|o ~ U/£q, where £q is a distance 
typically of the order of a few diameters (this point is 
confirmed in the full discussion to follow) . One thus gets 



T 



1/2 



p e. 



-m 



(38) 



which still depends on the ratio between shear stress a xy 
and pressure P. Combining this expression for To as a 
function of U and u xy /P and the formula for <J xy /P as 
a function of To in Eq. ([}?]) yields the final result 



'.ry 



P~ 



1-/3 



(39) 



where we have introduced the dimensionless constant, 

^ ~ \ Po d^e > ) • 

The cross-over between the two regimes occurs at the 
"critical" velocity U c defined as U c ~ {P jmpo) 1 ! 2 . 



3. Shear forces: quantitative 

These results obtained from the arguments of Sec. 
V.D.2 are in fact fully confirmed by a more detailed and 
careful analysis of the nonlinear set of equations for tem- 
perature and velocity Eq. (|l^) and Eq. ((2^). Here we 
present the full analysis. In order to simplify the discus- 
sion we assume that the granular material is semi-infinite, 
i.e. the cell width H is larger than the decay length for 
temperature 5 defined in Eq. (p3|). 

We proceed in two steps. First we obtain two closed 
implicit equations for the temperature at the boundary 
(To) and the shear stress a from an analysis of Eqs. ( |Tgj ) 
and (|20|). These equations are written in terms of the 
general form of the viscosity, Eq. ( |l4|) . Then in a second 
step, we discuss the different limiting behavior for the 
shear stress as a function of the shear velocity U ( "large" 
and "small" velocities). 

a. First condition denning To and a xv : We first 



rewrite the full temperature equation, Eq. (19), in terms 
of the expressions for the transport coefficients, Eqs. Jl4| ) 
and (III) : 



Ao 
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- tttT 1/2 = 0, 



(41) 



where we have introduced T s — 8 2 a 2 y pod A / (2X P) , a 
quantity with the dimension of a temperature; It can be 



rewritten more conveniently as T s = (P/po) (^Jr-) 2 
if we use the explicit expression of 6 which has been given 
in Eq. (p3[). It is important to note that the term fj(p(y)) 
is a function of T(y) only through the equation of state, 
Eq. (|9|). Let us introduce u = (T/m) 1 / 2 , which has the 
dimension of a velocity. A first integral of Eq. (fio|) can be 
obtained using a standard recipe of classical mechanics: 
we multiply Eq. ( fiO| ) by du/dy and integrate over u. We 
obtain: 



1 fdu 
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(42) 



Note that the density dependence of rj(p) has been rewrit- 
ten in terms of its dependence on u = (T/m) 1 / 2 thanks 
to the equation of state, Eq. (|l0|). The parameter E 
is a constant which in the mechanical analogy fixes the 
"energy" of the system. We now introduce the effective 
"potential" V(u) defined as 



V{u) 



du' 



T, 
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mS 2 r](pomu' 2 / P)u' 25 



^< (43) 



which allows us to rewrite Eq. (f42|) as 



+ V(u) = E. 



(44) 



The behavior of the potential V(u) can be obtained by 
analyzing the two limits u — > and u — > oo. The limit 
u ~ corresponds to small temperature, i.e. densities 
close to RCP. According to the expression of the viscosity 
close to RCP, Eq. @, the first term in the potential in 
Eq. (M3l) behaves in this limit like 
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(45) 



This term vanishes more rapidly than the last, quadratic 
term of V(u) (since the exponent verifies (3 > 1). There- 
fore close to u ~ 0, one has V(u) ~ — ^u 2 . 

In the limit of large u, i.e., high temperature and in- 
termediate density, the integral term in V(u) behaves as 
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which, when combined with the second term in Eq. (|4£ 
gives 
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Thus, depending on the sign of the first factor on the 
right had side of the equation, V{u) goes to plus or minus 
infinity when u — > oo. In Fig. O, we have arbitrarily 
chosen a positive sign to emphasize the different behavior 
of V(u) in the different limits, but the discussion does not 
depend on this particular choice. 




FIG. 11. Typical behavior of the potential V(u) as a func- 
tion of u. 



The "energy" E in Eq. (|42|) can be computed from the 
boundary conditions. As already mentioned, we assumed 
a semi-infinite system in this section to simplify the dis- 
cussion. At infinity, both the temperature and its gradi- 
ent are expected to vanish, and so does u : du/dy — > 
and u — > 0. Consequently the "energy" E is obtained to 
be zero : E — 0. This imposes the temperature at the 
boundary, T — mil 2 ,, which has to satisfy 



+ V(uo) = 0. 



(48) 



As we mentioned briefly in Sec. V.C, one expects in 
general that the boundary condition for the temperature 
relates the temperature gradient to the temperature at 
the boundary, in the form 



dT 

jQ = - K ^y 



RkT\ v= q, 



(49) 



y=o 



where Rk is a phenomenological parameter, usually de- 
noted as the Kapitza resistance. The ratio Ik — k/Rk 
has the dimension of a length. Following the same ar- 
guments which led us to the expression of the ther- 
mal conductivity, Eq. (|lg|), Rk is expected to be pro- 
portional to the fluid- wall collision frequency, i.e. to 
P/(poT 1/2 (y = 0)). The length £ K is thus expected to be 
independent of pressure, and only fixed by the "micro- 
scopic" quantities like the diameter of the particles and 
the roughness of the walls. Rewriting Eq. (M9) in terms 



of the field u, we find the condition for the temperature 
at the boundary 



^4 + v(u ) = o, 

with i = 21 k- Using Eq. @ for V(u), we obtain 
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i 2 , and T, 



where we recall that To = 

It is interesting to note that, since the second term in 
Eq. (|5l| ) is positive, the first term has to be negative, 
so that £ > S. If the opposite case, the only allowed 
solution to Eq. (|5l| ) is To = 0, i.e. no temperature is 
introduced in the system. This condition is in some sense 
expected : if I goes to zero, the boundary condition for 
the temperature at the wall, Eq. (f49|), imposes T = 0, 
so that no temperature can be introduced in the system. 
Physically I is related the amount of energy input at the 
boundary, while S has to do with the dissipation (see 
Eq. (^3[)), so the condition I > S can be understood as 
an energy balance condition : if not enough energy is 
introduced in the system at the boundary, the dissipation 
is too strong to allow flow. 



b. Second condition defining Tq and a xy 



A 



second condition relating a xy and To can be obtained by 
formally integrating the momentum equation, Eq. (EOT): 



V x (y) = U- 



r)(p,T) 



dy', 



(52) 



where we used the fact that a xy < to replace a xy by 
— \(J xy \ As y — > oo, the velocity goes to zero, so one gets 
the condition 



U 



dy'. 



(53) 



The explicit form of the y dependence of density and 
temperature is not known explicitly, but can be obtained 
formally from the first integral obtained previously, Eq. 
( [42| ) (with E = 0). The latter yields an implicit condition 
for the field u as 



dy = 



du 



(54) 



which allows us to rewrite the second condition, Eq. 
in terms of rj(p, T) only 
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(55) 



where the potential V(u) is given in Eq. (43). 

The two equations, Eqs. (pl|) and (55), in principle 



allow both To and a xy to be determined in terms of the 
pressure P and velocity U. 
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c. Limiting behaviors of the shear forces Two 

limits of the previous equations can be discussed for large 
and small shear velocity limit U (we shall specify below 
what we mean by "large" and "small" ) . 
(i) Large velocity limit 

Let us assume that the temperature at the boundary 
To satisfies To > P/po- This corresponds to the large 
velocity limit. In this case, the density close to the wall 
is intermediate and the shear viscosity takes its Enskog 
asymptotic expression, as quoted in Eqs. ( fuj) and (|l5|) : 
fjjpomu /P) ~ P/pomu 2 . 

In Eq. (pl|), corresponding to the first condition, the 
integral is dominated by its large u behavior, where the 
function fj behaves like fj(p mu 2 / P) ~ (P/mpo)u~ 2 . To 
a good approximation the first condition then takes the 
asymptotic form 
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which can be rearranged to give 
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Using the expression for T s , one gets 
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with p d = ^i] Q e /p Q d 3 (l - y)| 

Although the model we have developed is purely a hy- 
drodynamic model, this relationship is usually referred 
as to solid-like behavior. 

The temperature at the boundary can be obtained 
from the second condition, Eq. (|55|). The integral on 
the right hand side of Eq. (|55| ) is dominated by the be- 
havior near u ~ uq, with fj(pomu 2 / P) ~ P/pomu 2 and 
V(u) ~ V(u ) = -u\j2l 2 . One thus gets 
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which gives eventually 



T = mul = Q mil 2 , 



(59) 



(60) 



where the numerical prefactor ^ is defined as Cd — 

(2 Vo /(pdPod 2 £)) 1/2 . 

This shows that the large velocity regime which we 
defined as T > P/p ratifies U > (P / (p t mpo)) 1/2 ~ 
(Pd 3 /m). 

(ii) Small velocity limit 

In this case the viscosity is expected to behave "anoma- 
lously", with fj(p mu 2 /P) ~ rh(P / p mu 2 Y . In this 
limit, the first condition, Eq. (pl|), then takes the form 
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which can be reorganized to give 
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(62) 



This equation has to be completed by the second condi- 
tion, Eq. (|55|). In this equation, the potential term V(u) 



can be approximated by V(u) 



2S- 



■u , since the temper- 



ature T 



is small (see the discussion of the shape of 



the potential in the previous subsection). Together with 
the asymptotic expression of the viscosity in this limit, 
one can rewrite Eq. (pa) as 
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which yields 
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When combined with Eq. (|62j), one gets the following 
relationship between shear stress a xy and U : 



1-0 



•xy 



MwP{f/ 2 (^)}"- • 
The temperature at the boundary is then given as 



(66) 



where the two numerical prefactors pw and Cw are de- 
fined as p w = ((1 - ^)^p m) ( 2 ^mm { ^^-P)/P 

and C, w = ((2/3 - l)/(p wP od 2 6)) 1 /^- 1 \ The stress 
crosses over from its low-velocity power-law in velocity 
behavior to its high- velocity constant value at U « U c — 
P/(mp ). 

A crucial feature of the U dependence of shear stress 
in this regime is that it is velocity weakening, i.e., 
the shear stress decreases with increasing U for U < U c 
because (3 > 1. 



4- General remarks about shear forces at constant pressure 

Before ending this section, we would like to emphasize 
a few points : 

• we have identified two regimes in the shearing ve- 
locity dependence of the shear force. For small ve- 
locities, the shear force is velocity weakening, as ob- 
tained in Eq. (|65|). For large velocities, the shear 
force is found to be 11 solid-like" , i.e., independent of 
the shearing velocity, and proportional to pressure, 
as shown in Eq. (|58|). 
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• The cross-over between one regime and the other 
is found to occur at a critical velocity, U c = 
(Pd 3 /m)^ 2 . 

• In the velocity weakening regime, i.e. for veloci- 
ties smaller than U e , the steady sliding regime is 
expected to be unstable. Indeed, as da xy /dU < 0, 
the steady sliding situation is unstable to small fluc- 
tuations (see Refs. [p2p^| for a full discussion). As 
a result, one might expect a stick slip behavior of 
the system for U < U c . 

• Finally, it is interesting to note that both the large 
and small U limits of the shear rates, Eqs. ( J5q ) 
and (|65|), are in agreement with our previous di- 
mensional analysis of the dynamics in the constant 
volume case (see Sec. V.D.I). Since on dimensional 
grounds the pressure should be proportional to U 2 
when the volume is fixed, one finds that both ex- 
pressions reduce to the Bagnold scaling a xy oc U 2 . 
We emphasize that this "global" scaling holds even 
if the local Bagnold relationship a xy oc j 2 is not ex- 
pected to hold in the small velocity regime (where 
the viscosity exhibits an "anomalous" behavior). 

VI. CONCLUSIONS 

In this paper, we have investigated both experimen- 
tally and theoretically the shear dynamics of granular 
flow. The granular material was sheared in a Couette 
cell, instrumented to allow measurements of both shear 
forces and "microscopic" dynamics of individual parti- 
cles. The latter measurements were performed by track- 
ing the instantaneous positions of particles on the upper 
surface using a fast camera and an imaging procedure. A 
variable upward air flow through the granular material al- 
lows the internal pressure to be adjusted. Experimentally 
the following results were found: (i) the flow is restricted 
to a small region close to the moving boundary; (ii) the 
normalized velocity profile is independent of the shearing 
velocity, pressure (i.e. airflow), and of the type of motion 
of the moving cylinder (stick-slip or continuous sliding); 
(iii) the shear force acting on the moving cylinder is in- 
dependent of the shearing velocity and proportional to 
pressure, as shown in the planar geometry in j33|); (iv) 
there is a close connection between local fluctuations and 
mean flow, as manifested in the power law relationship 
measured experimentally between these two local quan- 
tities (Fig. 10); (v) the RMS velocity profile decays more 
slowly as a function of the radial coordinate than does 
the mean velocity profile. We find that these general fea- 
tures occur for several different types of particles, though 
there are quantitative differences. 

On the basis of these observations, we have proposed 
a "hydrodynamic" model for the granular flow, in which 
the granular material is assumed to behave like a locally 



Newtonian fluid. In contrast to "classical" fluids how- 
ever, the temperature and density of the material are 
not constant over the shear-cell. The temperature, de- 
fined in terms of the fluctuations of the velocity, is cre- 
ated at the moving boundary and propagated through 
the material. The temperature profile thus results from 
a balance between heat flux and local energy loss due 
to the inelasticity of collisions. As a consequence, the 
density and temperature dependences of the transport 
coefficients play a crucial role. 

In the high density regime under consideration, sim- 
ple, asymptotic expressions for the transport coefficients 
can be obtained within the Enskog approximation. The 
latter is however expected to be invalid for the density 
dependence of the viscosity at very high densities (close 
to Random Close Packing, RCP) where collective rear- 
rangements comes into play. This led us to assume a 
stronger divergence of the viscosity as a function of the 
density close to RCP. In a manner analogous to what is 
usually proposed in supercooled liquids, we have assumed 
an algebraic divergence of the viscosity near the RCP 
density. Using the equations of transport of momentum 
and heat, we then compute the mean and RMS velocity 
profiles. Those are found to be in very good agreement 
with the experimental results. Moreover, the scaling law 
relationship between the mean and RMS velocity pro- 
files that is found experimentally is also predicted from 
the model. 

Finally, the velocity dependence of the shear force is 
determined. Two regimes are predicted, depending on 
whether the velocity is larger or smaller than a critical 
velocity U c = (Pd 3 /m) 1 / 2 . For large velocities a solid like 
friction force is predicted, i.e. independent of the shear- 
ing velocity, in agreement with the experimental obser- 
vation. For small velocities, a velocity weakening regime 
is predicted, consistent with the occurrence of stick-slip 
motion in this case. We emphasize that these results are 
found within the hydrodynamic model, even if no solid 
friction force is assumed to hold between the grains. 

In our work, we have assumed that the flow on the 
upper surface is close to that within the interior of the 
Couette cell, as found in earlier work ||. However, 
since the particles on the surface are less constrained, 
their fluctuations could be somewhat different from those 
in the interior, even when we apply a downward pres- 
sure through an airflow. While MRI and X-ray mea- 
surements in the interior || could not resolve fluctua- 
tions, recent unpublished measurements of fluctuations 
on the more constrained bottom surface layer by the 
Chicago group reveal fluctuation profiles similar to 
ours: Fluctuations decay roughly exponentially with a 
characteristic length longer than that for the mean veloc- 
ity. The velocity profiles found in these two experiments 
differ somewhat; however differences in geometry (inner 
and outer radii), material properties, and measurement 
methodology make a direct comparison between the ex- 
periments difficult. 

The main point of the present paper is the agreement 
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between our measured velocity and temperature profiles 
(and friction force) and the predictions of a Newtonian 
hydrodynamic model with a strong density dependence 
of the viscosity taken from models of the glass transition. 

Many questions are still open. One important problem, 
which is not discussed here, is the temporal response of 
the system when a velocity step is imposed. Such an 
experiment would provide information about the trans- 
port mechanism within the granular material. The re- 
sponse can be predicted in principle from the hydrody- 
namic model and could provide an independent test of 
its validity. 
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